Genomic analysis reveals deep population divergence in the water snake Trimerodytes percarinatus (Serpentes, Natricidae)

Abstract Although several phylogeographic studies of Asian snakes have been conducted, most have focused on pitvipers, with non‐venomous snakes, such as colubrids or natricids, remaining poorly studied. The Chinese keelback water snake (Trimerodytes percarinatus Boulenger) is a widespread, semiaquatic, non‐venomous species occurring in China and southeastern Asia. Based on mitochondrial DNA (mtDNA) and single nucleotide polymorphism (SNP) data, we explored the population genetic structure, genetic diversity, and evolutionary history of this species. MtDNA‐based phylogenetic analysis showed that T. percarinatus was composed of five highly supported and geographically structured lineages. SNP‐based phylogenetic analysis, principal component analysis, and population structure analysis consistently revealed four distinct, geographically non‐overlapping lineages, which was different from the mtDNA‐based analysis in topology. Estimation of divergence dates and ancestral area of origin suggest that T. percarinatus originated ~12.68 million years ago (95% highest posterior density: 10.36–15.96 Mya) in a region covering southwestern China and Vietnam. Intraspecific divergence may have been triggered by the Qinghai‐Xizang Plateau uplift. Population demographics and ecological niche modeling indicated that the effective population size fluctuated during 0.5 Mya and 0.002 Mya. Based on the data collected here, we also comment on the intraspecific taxonomy of T. percarinatus and question the validity of the subspecies T. p. suriki.

Intraspecific divergence may have been triggered by the Qinghai-Xizang Plateau uplift.Population demographics and ecological niche modeling indicated that the effective population size fluctuated during 0.5 Mya and 0.002 Mya.Based on the data collected here, we also comment on the intraspecific taxonomy of T. percarinatus and question the validity of the subspecies T. p. suriki.

| INTRODUC TI ON
The patterns of genetic diversity, population structure, and evolutionary history of organisms have been greatly affected by various biological factors and abiotic factors.Of those, geography and climatic oscillations have been previously shown to generate population structure, drive the process of speciation, and influence population size contraction and expansion (Burbrink et al., 2016(Burbrink et al., , 2022;;Hewitt, 2004;Myers et al., 2017;Qu & Lei, 2009;Yang et al., 2022).
Complex geological events can result in physical barriers, such as mountains, rivers, and marine incursions, which can prevent gene flow between closely related species or populations.Furthermore, unfavorable climate changes should cause populations to contract, on the contrary, a favorable climate may lead to population expansion (Berv et al., 2021;Fijarczyk et al., 2011;Li et al., 2022;Myers et al., 2019Myers et al., , 2020;;Pyron & Burbrink, 2009;Ursenbacher et al., 2008;Yang et al., 2022;Zhao et al., 2022).However, the influence of both factors on these processes differ in the timing and effect among different regions (such as Asia and North America) and among groups (such as birds and reptiles).
The Chinese keelback water snake, Trimerodytes percarinatus, is a semiaquatic, nonvenomous colubrid with a wide geographic distribution in China and Southeast Asia (Das, 2010;Uetz et al., 2023;Zhao, 2006).It is primarily nocturnal and can be found in various waterbodies, such as streams, rivers, and ponds, at elevation ranging from 100 to 2000 meters above sea level (Zhao, 2006).Currently, two subspecies are generally recognized, Trimerodytes percarinatus percarinatus (Boulenger, 1899), based on type specimens from Fujian, China, and Trimerodytes percarinatus suriki (Maki, 1931), based on specimens from Taiwan, China (Zhao et al., 1998).The two subspecies differ from each other only in terms of the number of ventral scales (Zhao et al., 1998).However, Zhao (2006) questioned the validity of T. p. suriki and suggested that further research was needed to evaluate its taxonomic status.
The widespread geographic distribution and specialized ecology make T. percarinatus a good model for evaluating the impact of geographical features and climatic cycles on the demographic history of species in East and Southeast Asia.Herein, based on extensive geographic sampling, integration of mitochondrial DNA (mtDNA) and genome-wide single nucleotide polymorphism (SNP) markers, we explored the population structure and evolutionary history of this species.We aimed to (1) determine the population structure and evaluate genetic diversity and (2) establish the timing of species diversification and evolutionary dynamics within the species.

| Taxon sampling and DNA extraction
Samples were obtained via fieldwork or provided from colleagues and natural history museums.In total, 180 samples of T. percarinatus from 71 localities covering almost its entire geographic range were collected, sequenced, and analyzed (Figure 1; Table S1).In addition, five samples from known congeners T. annularis, T. balteatus, T. aequifasciatus, T. yunnanensis, and T. yapingi were included, with T. yapingi chosen as the outgroup based on a previous phylogenetic study (Guo et al., 2020).Total genomic DNA was isolated from liver, muscle, or skin tissue preserved in 85% ethanol using an Animal Genomic DNA Purification Kit (Tiangen Bio-Tech Co., Ltd., Beijing, China).

| Mitochondrial DNA sequencing and SNP generation
Two mitochondrial gene fragments, cytochrome b (cyt.b) and NADH dehydrogenase subunit 2 (ND2), were amplified with the primers L14910/H16064 (Burbrink et al., 2000) and hxND2F/hxND2R (Reynolds, 2011), respectively.The polymerase chain reaction (PCR) parameters were identical to those described in the corresponding references for each primer pair.The PCR products were purified, The raw RAD-seq Illumina reads were subjected to quality control and filtered using the Novogene QC program (parameter: -r 150 -N 0.1 -q 33 -L 5 -p 0.5) to obtain reliable reads and avoid those with artificial biases.The following types of reads were removed during quality control: (i) those containing ≥10% unidentified nucleotides; (ii) those with more than 10 nucleotides aligned to the adaptor or mismatches exceeding 10%; (iii) those with over 50% of read bases with Phred quality scores lower than 5; and (iv) putative PCR duplicates generated during library construction.The RAD-seq data were assembled using SOAPdenovo2 (Luo et al., 2012) with the parameters (max_rd_len = 150, avg_ins = 200, reverse_seq = 0, asm_flags = 3, rd_len_cutoff = 150, rank = 1-K 47) to generate contigs.The contigs were aligned to the pseudo-reference sequence by using BWA (Li & Durbin, 2009) with the command "mem -t 4 -k 32 -M" to generate SAM files, which were imported into SAMtools (Li et al., 2009) for sorting and merging, and then converted into BAM files.Simultaneously, both sequencing depth and coverage were calculated according to the BWA alignment results.Finally, BAM files were imported into SAMtools (Li et al., 2009) for SNP detection.
The minimum percentage of individuals in a population required to retain a locus for that population (-r) was set to 80%.All SNPs with a minor allele frequency (maf) < 0.05 were excluded from the genotype dataset of all samples.

| Phylogenetic analyses
All mtDNA sequences were assembled and edited manually using SeqMan in Lasergene v15.1 (DNASTAR, USA).Sequences were aligned using the MUSCLE algorithm implemented in MEGA v7 (Kumar et al., 2016), and both segments were concatenated in the software MEGA v7 (Kumar et al., 2016).Sequences were converted into amino acid residues to check for stop codons and ensure an open reading frame.PartitionFinder v2.1 (Lanfear et al., 2017) was used to identify the best partitioning schemes and model of sequence evolution for each partition based on the Akaike information criterion (AIC).Bayesian inference (BI) and Maximum likelihood (ML) methods were used to infer phylogenetic relationships.
Phylogenetic analyses were conducted for the SNP dataset using a ML approach with the GTRGAMMA model and outgroup T. yapingi.The ML analysis was performed in RAxML v8.2.11 (Stamatakis, 2006), with branch support assessed by 100 rapid bootstrap replicates.
To assess the number of genetic groups within T. percarinatus, we used multivariate and model-based methods for the SNP data.
First, a principal component analysis (PCA) was conducted in Plink v1.9 (Chang et al., 2015), and the first three principal components were plotted.Second, population structure was inferred assuming different numbers of clusters (K) from 2 to 8 using Admixture v1.3.0 (Alexander et al., 2009).We used 10-fold cross-validation (CV) to compare different numbers of clusters, with the lowest CV value indicating the most likely number of clusters.

| Gene flow analyses
To investigate gene flow among the SNP-based lineages, two different methods were used.First, the D and f4 statistics (from the ABBA-BABA test; Green et al., 2010) were used to estimate intraspecific gene flow between all of the T. percarinatus lineages based on SNP data using "qpDstat" in AdmixTools (Patterson et al., 2012), T. yapingi was also included as the outgroup.The statistical significance of D and f4 statistics was assessed through a block-jackknifing approach as implemented in AdmixTools, in which the absolute value of a Zscore > 3 was considered significant (Patterson et al., 2012).Second, TreeMix v 1.13 (Pickrell & Pritchard, 2012), which applies an ML method based on a Gaussian model of allele frequency change, was conducted using 5 different migration parameter settings, ranging from zero to 4 migration events.For each migration setting, we conducted 1000 bootstrap replicate analyses and identified the highest likelihood run from this set of 1000 to characterize evidence for the best-fit migration model.

| Divergence date and ancestral area estimation
BEAST v 2.6.6 (Bouckaert et al., 2014) was used to estimate the date of origin of T. percarinatus as well as each lineage based on all mtDNA sequences using a Yule model and a relaxed uncorrelated lognormal clock (Drummond et al., 2006).The following three calibration constraints were used: (1) The most recent common ancestor (MRCA) of all Natricinae was constrained with a lognormal mean of 30 Mya and standard deviation (SD) of 0.115, yielding a 95% confidence interval (CI) from 23 to 36 Mya; (2) MRCA of Natrix was constrained with a lognormal mean of 22 Mya and SD of 0.15, yielding a 95% CI of 16-28 Mya; (3) MRCA of the genus Thamnophis was constrained with a lognormal mean of 16 Mya and SD of 0.15, yielding a 95% CI interval of 12-20 Mya (Guo et al., 2012).For these calibrations, several additional sequences were retrieved from GenBank, with Acrochordus granulatus chosen as the outgroup (Table S1).Two independent searches were conducted based on 2 × 10 7 generations and sampling every 1000 iterations, with the first 25% of samples discarded as burn-in.Both analyses produced similar results and yielded ESS values >200 for the posterior probability distribution of all parameters (Drummond et al., 2006), as calculated in Tracer v 1.7 (Rambaut et al., 2018).
To infer the biogeographic history of T. percarinatus, we estimated ancestral area using dispersal vicariance analysis (S-DIVA) and dispersal extinction cladogenesis (DEC) using the program Reconstruct Ancestral States in Phylogenies (RASP; Yu et al., 2014).
Analyses were conducted for one million cycles, with 25% of the initial run discarded as burn-in.The reconstructed biogeographic region for each node was summarized and plotted as a pie chart.

| Demographic history
To determine changes in population size over time, we examined the historical demographics within T. percarinatus lineages using Multiple Sequentially Markovian Coalescent (MSMC2, Schiffels & Wang, 2020) based on the SNP dataset.First, the program bcftools (Li, 2011) was used to filter out the sites with missing genotype ". /." from each lineage of the SNP-based phylogenetic tree, and two samples were randomly selected for each lineage.Then, the MSMC2 analysis was executed following Dong et al. (2023) with the parameters as follows: (1) the generation time (g) was set to 3 years; (2) the generation mutation rate (μ) was 4.71 × 10 −9 per site per generation (Peng et al., 2020).

| Ecological niche modeling
To evaluate the influences of glacial cycles on the geographic distribution of T. percarinatus, we predicted suitable habitat during the Last Glacial Maximum, Mid-Holocene, and at present using ecological niche modeling (ENM).Geographic coordinates were obtained from a variety of sources, including museum records of 535 T. percarinatus specimens from Vertnet (www.vertn et.org), published accounts (e.g., Yang & Rao, 2008;Zhao et al., 1998), and field surveys.After filtering duplicate locations that were within 5 km 2 of one another, a total of 233 geographic coordinates were retained for T. percarinatus.The 19 bioclimatic variables for the present , Mid-Holocene (~0.006Mya), and the Last Glacial Maximum (LGM, 0.022-0.011Mya) were downloaded from WorldClim v1.4 (http:// www/ world clim.org) with a resolution of 2.5 arc minutes (Hijmans et al., 2005).We used a Pearson's correlation coefficient threshold of 0.85 to identify correlated variables.We retained 11 climatic variables (Table S2) after removing 8 climatic variables with high correlations.The ENM was conducted using Maxent v 3.4.4(Phillips et al., 2006).We used 75% of coordinates as training data and 25% as testing data.Simulations were conducted for 10 replicates, and the average result was selected as the best model output.The area under the receiver operating characteristic curve (AUC) was used to evaluate simulation accuracy, with AUC values >0.9 indicating good performance.We used Arcgis v10.4 (ESRI) to visualize the ENMs at the LGM, Mid-Holocene, and present day.
After assembling and filtering these genomic data, a total of 393,329 SNPs were obtained (mean = 6051 SNPs per sample) from all samples.

| Phylogenetic analyses
The mtDNA-based BI and ML gene trees showed similar topologies, with slight disagreements at some shallow nodes (Figure 2).Analyses consistently indicated that all T. percarinatus samples formed a wellsupported monophyletic group (1.00 PP and 100% UFB) with five major lineages, each receiving high support values (A-E; Figures 1   and 2).Lineage A, which was sister to all other lineages, consisted The SNP-based phylogenetic relationships were generally consistent with those estimated from mtDNA, but differed in a few specific ways.The ML analyses indicated that all T. percarinatus samples formed a highly supported monophyletic group (BS 100%) with only four distinct and well-supported lineages (I-IV; Figure 3).The samples from mtDNA lineages A, C, and D formed distinct and highly supported lineages labeled as I, III, and IV with the genomic data.
While those samples from lineages B and E formed a monophyletic lineage II.It is likely that these lineages are parapatric.The interlineage relationships were different in both analyses, but well resolved in the genomic data tree.

| Genetic diversity and population structure
Based on the SNP data, Ho ranged from 0.1103 (lineage IV) to 0.1587 (lineage I) with a mean of 0.1433 per lineage; He ranged from 0.1458 (lineage IV) to 0.2247 (lineage II), with an average of 0.1815.The estimated Fis ranged from 0.0289 (lineage I) to 0.3337 (lineage II).
The highest nucleotide diversity was found in lineage II (0.0770) and the lowest in lineage IV (0.0041) (Table 1).The uncorrected interlineage genetic distances (p-distances) ranged from 2.88% (III and IV) to 26.81% (I and II); Pairwise Fst values were all statistically significant, ranging from 0.1633 (II and III) to 0.3428 (I and IV) (Table 2).
The mean interlineage uncorrected p distances ranged from 3.39% to 6.02% based on cyt.b and from 2.77% to 4.86% based on ND2 (Table S3).The overall Hd was 0.96, π was 0.029-0.031(Table S4).4).Admixture analyses found that K-values ranging from 3 to 5 fit the data (Figure S1).When K = 4, the structure corresponded to the topology of the concatenated SNP gene tree (I-IV; Figures 3 and 5).The admixture analyses also demonstrated high levels of admixed ancestries for several sampled individuals.

| Gene flow analyses
The results from the D and f4 statistics consistently revealed that all tests significantly deviated from neutrality (|Z-score| > 3), suggesting the presence of gene flow between lineages I and III, as well as lineages II and IV (Table S5).In addition, comparing the likelihood  support of the 5 different TreeMix runs suggested that there was only 1 instance of migration in the population tree model, which occurred between lineages I and III (Figure 6).

| Historical demography
MSMC2 analyses revealed continuous declines of effective population size (Ne) from 0.5 to 0.05 Mya for each analyzed lineage of T. percarinatus (Figure 8).Subsequently, Ne of the four lineages remained relatively stable from about 0.05 to 0.005 Mya.Then, all four lineages experienced a drastic rise in effective population size around 0.006 to 0.002 Mya, followed by a second sharp decline (Figure 8).

| Ecological niche modeling
The AUC value approached 1 (≥0.98),indicating a good performance of the predictive models.The ENM results revealed that during the LGM, the Sichuan Basin, southern Guangdong and Guangxi, the south and east of Taiwan in China, northern Vietnam and northwestern Myanmar contained suitable environments and may have served as potential refugia (Figure 9A).During the Mid-Holocene, suitable habitats extended southward to southern Vietnam and northward to the Yangtze River region in China (Figure 9B).The overall predicted distribution pattern was largely consistent with the current distribution of T. percarinatus (Figure 9C).

| Phylogeographic patterns across East Asia
Snakes are generally susceptible to climate fluctuation and associated habitat changes over time.Although snake phylogeographic studies in China lag far behind those in North American and Europe, a few studies have revealed phylogeographic patterns across China (Guo et al., 2011(Guo et al., , 2016(Guo et al., , 2019;;Huang et al., 2007;Zhu et al., 2016).However, these studies are mainly based on venomous snakes, particularly pitvipers.Our present study is the first work using a nonvenomous snake species with range-wide sampling and genome-based data.The results presented here demonstrate conclusively that T. percarinatus is composed of multiple ancient geographic lineages distributed across China and Vietnam.and Hengduan Mountains Region were the center of origin (Guo et al., 2011(Guo et al., , 2016(Guo et al., , 2020;;Zhu et al., 2016).Our results for T. percarinatus are consistent with both this hypothesis and "out of Qinghai-Xizang Plateau" and "out of Hengduan Mountains" (Guo et al., 2019;Päckert et al., 2020).The estimated origin date of this species is 12.68 Mya in the mid-late Miocene, which is much older than some codistributed pitvipers (Guo et al., 2011(Guo et al., , 2016)), but similar to other colubrids (Guicking et al., 2006).Importantly, this date falls within the estimated time frame of the QXP uplift during the Miocene (25-10 Mya; Sun, 1997).The uplift of the QXP resulted in environmental transitions and historical climate change, and therefore may have facilitated T. percarinatus population divergence as seen in other taxa (Che et al., 2010;Fu et al., 2005;Guo et al., 2020;Zhu et al., 2022).F I G U R E 7 Divergence date estimation of Trimerodytes percarinatus using BEAST v2.6.6 with mtDNA loci.All genetic lineages are color-coded as Figure 1.
In East Asian snakes, phylogeographic studies have typically revealed the presence of both longitudinal and latitudinal diversification (e.g., Guo et al., 2011Guo et al., , 2020;;Huang et al., 2007;Lin et al., 2014).Guo et al. (2020) suggested that longitudinal diversification occurred more commonly than latitudinal divergence across this region.Phylogenetic analyses based on the SNP dataset show that both diversification patterns were detected within T. percarinatus.For example, latitudinal diversification had occurred between lineages I + II and III + IV (Figure 3).The Yunnan-Guizhou Plateau, which is located between I + II and III + IV, may be responsible for this diversification.Lee (1939) and Jiang and Wu (1993)  Mya (Shi et al., 2006;Zhao et al., 2007;Zhu, 2016).However, the estimated population divergence date indicates that the Hainan lineage of T. percarinatus originated ~6 Mya, suggesting that diversification predates the isolation of this island.Similar cases have been reported in several species of Asian pitvipers (Guo et al., 2016(Guo et al., , 2019)).

| Historical demography of T. percarinatus
Historical demography can be affected by various factors, including biological, such as life history, competitors, predation, parasitism, disease, human activities, and abiotic, such as geological processes, and glacial and interglacial periods.In Europe and North America, many organisms experienced changes in population size in response to the LGM (Guiher & Burbrink, 2008;Hewitt, 2004;Hull & Girman, 2005;Schield et al., 2018); however, this trend was not observed in numerous taxa distributed across East Asia (Gao et al., 2012;Guo et al., 2016Guo et al., , 2019;;Huang et al., 2013;Zhou et al., 2013).Here, analyses of demographic history indicated that the effective population size experienced great changes in recent date.In the Middle Pleistocene (about 0.5-0.05Mya), the formation of the Qinghai-Xizang Plateau resulted in inland arid zone in Asia (Zhang, 2011).This unsuitable habitat has led to population declines in T. percarinatus lineages.In the subsequent LGM of the late Pleistocene (about 0.05Mya-0.01Mya),the population sizes remained stable, indicating the effect of the LGM on the effective population size of T. percarinatus was not significant, which may be due to the refugia provided in the south region of the Yangtze River.ENMs confirmed the existence of several refugias in southwest China, southern China, and Southeast Asia (Figure 9A).et al., 2013;Wu et al., 2023).

| Intraspecific divergence within T. percarinatus and taxonomic implications
A high degree of intraspecific genetic diversity can be indicative of a species' ability to respond and adapt to changing habitats and environmental stressors.In this study, genomic sequence data demonstrated high genetic diversity, and the presence of several geographically structured lineages within T. percarinatus.The high genetic diversity and pronounced population structure in T. percarinatus may be attributed to several different factors.As a semiaquatic species, which prefers small rivers, streams, and ponds, the dispersal ability of T. percarinatus is contingent upon drainage systems and river networks.Therefore, the lack of these aquatic features in terrestrial ecosystems may act as a physical barrier to gene flow between populations, thereby promoting genetic differentiation (e.g., Yan et al., 2013;Zhou et al., 2013).Population structure analyses (Figure 5) indicated that there was no significant gene flow among lineages.Second, T. percarinatus typically inhabits regions at elevations below 1000 m (rarely exceeding 2000 m), and therefore mountain ranges may also prevent dispersal.
The population that occurs on Taiwan, long recognized as the subspecies T. p. suriki (Maki), was represented in this study by a single sample (GP 5122).Despite the lack of SNP data generated for this sample, the mtDNA-based phylogeny indicated that it was nested within the clade of all samples collected from the Chinese mainland, resulting in the nonmonophyly of the nominal subspecies.
Furthermore, the Taiwanese T. p. suriki ventral scale counts (range: 142-153) are not significantly divergent from its closest relatives T. p. percarinatus (range: 131-160) in China (Zhao et al., 1998).Thus, based on this evidence, we question the validity of the subspecies T. p. suriki, and suggest that further investigation is needed.

| 3 of 15 LYU
et al. with the double-stranded products bidirectionally sequenced at a commercial sequencing company (Sangon, Chengdu, China).Based on the mtDNA phylogenetic analysis (see below), 61 individuals representing all mtDNA lineages were selected to generate restriction site-associated DNA sequencing (RAD-seq) libraries following the protocols of Peterson et al. (2012) and Schield et al. (2015).In brief, DNA was digested with the restriction enzyme EcoRI, and unique barcode adapters were ligated onto each individual sample in the library.Sequencing of 200-400 bp fragments was performed on a single lane of an Illumina HiSeq sequencer using a paired-end protocol by Novogene Ltd. (Tianjin, China).

F
Map showing localities of Trimerodytes percarinatus sampled in this study.Symbols indicate different mtDNA lineages and colors represent SNP-based lineages.Square: A, Triangle: B, Circle: C, Pentagram: D, Diamond: E. Blue: I, Red: II, Yellow: III, Green: IV.
of the specimens from Hainan Island exclusively.Lineage B was composed of specimens from southwestern Guangxi in China and northeastern Vietnam.Lineage C contained samples from eastern and central China, including Taiwan, Hong Kong, Guangdong, Fujian, Zhejiang, Jiangxi, Anhui, southern Henan, eastern Hunan, and northeastern Guangxi.One specimen from Taiwan was sister to all other samples in lineage C. Lineage D included specimens from southwestern China, including Sichuan, Chongqing, Guizhou, western Hubei, western Hunan, and northwestern Guangxi.Lineage E, sister to lineage D, contained samples from Yunnan in China and the rest of Vietnam, excluding northeastern Vietnam.Although lineages B-E formed highly supported groups, the relationships among B, C, and D + E were unresolved (Figures 1 and 2).
Results of the PCA and admixture analysis based on SNP data supported four distinct lineages within T. percarinatus showing strong geographic structure.In the PCA, the first, second, and third principal components (variances explained: 14.18%, 11.41% F I G U R E 2 Bayesian 50% majority-rule consensus tree of Trimerodytes percarinatus inferred from the mtDNA dataset.Posterior probabilities from BI (>50%) and bootstrap support values from ML analyses (>50) are adjacent to major nodes.Branch support indices are not given for most shallow nodes to preserve clarity.Square: A, Triangle: B, Circle: C, Pentagram: D, Diamond: E. Blue: I, Red: II, Yellow: III, Green: IV. and 9.14%, respectively) separated this species into four different clusters (Figure

F
I G U R E 3 Maximum likelihood phylogeny inferred from genome-wide SNP data.Numbers above branches indicate bootstrap values for each node.Square: A, Triangle: B, Circle: C, Pentagram: D, Diamond: E. Blue: I, Red: II, Yellow: III, Green: IV.TA B L E 1 Genetic diversity of SNP data for each lineage.

F
I G U R E 4 PCA showing population genetic structure and clustering among Trimerodytes percarinatus individuals based on SNP data.All genetic lineages are color-coded as Figure 1.LYU et al. Results from ancestral area estimation indicate that T. percarinatus originated in an area located in southwestern China and Vietnam.Previous studies based on snakes at different taxonomic levels consistently suggested that Qinghai-Xizang Plateau (QXP)

F
I G U R E 5 Population structure and admixture proportions of Trimerodytes percarinatus based on SNP data.All genetic lineages are colorcoded as Figure 1.F I G U R E 6 Treemix admixture graph with one migration event (m) from SNP data.The migration weight and the directionality of gene flow are indicated by the colored arrows.The graded bar on the right approximates color to values of migration weight as the fraction of ancestry derived from the migration source, so that the warmer color (dark red) indicates an ancestry of 50%.The scale below measures the genetic drift from the ancestral to the extant populations.
suggested that the Asian continent was composed of three geomorphologic steps, gradually decreasing in elevation from west to east.The longitudinal diversification pattern observed here between lineages III and IV is consistent with the boundary between the second and third geomorphologic steps.Similar patterns have also been uncovered in three Asian pitvipers (e.g., Deinagkistrodon acutus: Huang et al., 2007; Protobothrops mucrosquamatus: Guo et al., 2019; Viridovipera stejnegeri: Guo et al., 2016).These similar results suggest that differences in habitat between the second and the third geomorphologic steps may lead to population divergence within codistributed species.Of course, more phylogeographic studies on Chinese snakes with similar geographic distributions will provide additional insights into the drivers of lineage formation.The geographic isolation of Hainan Island promoted lineage divergence between the local T. percarinatus population and those dispersed across mainland Southeast Asia.Hainan Island became geographically isolated from the Asian continent approximately 2.5 Both SNP-based phylogenetic analysis and admixture analysis indicate that the Hainan individuals are most closely related to the individuals from southern Guangxi, China and Vietnam, suggesting that the Hainan individuals may have been colonized from these regions, consistent with previous studies on other animal and plant groups(Chen et al., 2015;Lin et al., 2010;Zhu, 2016).Based on their work on snakes and previous studies on other organisms inhabiting in Hainan Island, Guo et al. (2016) proposed that there were two dispersal events contributing to the Hainan fauna and flora (Two Dispersals Hypothesis): one is from South China via the Qiongzhou Strait and the other from Vietnam via the Gulf of F I G U R E 9 Inferred present and historical suitable habitats based on the ENMs of Trimerodytes percarinatus localities.Suitable habitats for Trimerodytes percarinatus during Last Glacial Maximum (a), Mid-Holocene (b), and present (c) estimated using ENM.Warmer colors indicate regions with a higher logistic threshold probability of suitable habitat, while dark colors represent areas with lower probability of Trimerodytes percarinatus presence.F I G U R E 8 Demographic history of Trimerodytes percarinatus from SNP data.Time (in years) is shown along the x axis and effective population size is shown along the y axis.The generation time (g) is set to 3 years, and the generation mutation rate (μ) is 4.71 × 10 −9 per site per generation.All genetic lineages are colorcoded as Figure 1.LYU et al.Tonkin.The case of Hainan individuals colonizing from Vietnam via the Gulf of Tonkin further confirms the hypotheses on Hainan Island fauna proposed by Guo et al. (2016).